# CLASS 1 - SCRIPT #2
# (C) Fernando Ascensão 2025

library(geodata)

mywd <- "C:/Users/fjascensao/OneDrive - Universidade de Lisboa/Aulas/MIBC/Classes/Class_1_20260302"  # <- change to your root!

setwd(mywd)
dir()

mystudyarea <- st_read("Administrative/Iberia.shp")


# Worldclim 
?worldclim_global
mypath = file.path(mywd, "Env_var")


# includes Portugal and Spain / will create a new 'climate' folder
worldclim_country(country="ES", var="bio", path=mypath, version="2.1")

setwd(file.path(mypath, "climate", "wc2.1_country"))
dir()

bioclim_es <- rast("ESP_wc2.1_30s_bio.tif")

plot(bioclim_es)


# They are coded as follows:
#   
# BIO1 = Annual Mean Temperature
# BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))
# BIO3 = Isothermality (BIO2/BIO7) (×100)
# BIO4 = Temperature Seasonality (standard deviation ×100)
# BIO5 = Max Temperature of Warmest Month
# BIO6 = Min Temperature of Coldest Month
# BIO7 = Temperature Annual Range (BIO5-BIO6)
# BIO8 = Mean Temperature of Wettest Quarter
# BIO9 = Mean Temperature of Driest Quarter
# BIO10 = Mean Temperature of Warmest Quarter
# BIO11 = Mean Temperature of Coldest Quarter

# BIO12 = Annual Precipitation 
# BIO13 = Precipitation of Wettest Month
# BIO14 = Precipitation of Driest Month
# BIO15 = Precipitation Seasonality (Coefficient of Variation)
# BIO16 = Precipitation of Wettest Quarter
# BIO17 = Precipitation of Driest Quarter
# BIO18 = Precipitation of Warmest Quarter
# BIO19 = Precipitation of Coldest Quarter



# crop to study area

mystudyarea4326 <- st_transform(mystudyarea, 4326)


bioclim_iberia <- terra::crop(bioclim_es, mystudyarea4326, mask=TRUE)

plot(bioclim_iberia) # remove 8 and 9
plot(bioclim_iberia[[-c(8,9)]])


# change resolution and project to 3035
my_bioclim <- project(bioclim_iberia[[-c(8,9)]], crs("EPSG:3035"), res=1000)
names(my_bioclim) <- gsub("wc2.1_30s_", "", names(my_bioclim))


dir.create(file.path(mywd, "Env_var","Bioclim"))
setwd(file.path(mywd, "Env_var","Bioclim"))  #  <--- attention where you store this!!!!!!

writeRaster(my_bioclim, "my_bioclim.tif", overwrite=T)
writeRaster(my_bioclim, filename=paste0(names(my_bioclim), ".asc"), overwrite=T, NAflag=-9999)



# future climate ----------------------------------------------------------

library(pastclim)


list_available_datasets()[grepl("WorldClim_2.1",list_available_datasets())]


get_vars_for_dataset(dataset = "WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m")


set_data_path(path_to_nc =file.path(mywd, "Env_var", "Future_climate"))


# # Future projections are based on the models in CMIP6, downscaled and de-biased using WorldClim 2.1 for the present as a baseline.
# Monthly values of minimum temperature, maximum temperature, and precipitation, as well as 19 bioclimatic 
# variables were processed for 23 global climate models (GCMs), 
# and for four Shared Socio-economic Pathways (SSPs): 126, 245, 370 and 585. 
# Model and SSP can be chosen by changing the ending of the dataset name WorldClim_2.1_GCM_SSP_RESm.

# # A complete list of available combinations can be obtained with:
# 
# # list_available_datasets()[grepl("WorldClim_2.1",list_available_datasets())]
# # So, if we are interested in the the HadGEM3-GC31-LL model, with ssp set to 245 and at 10 arc-minutes, we can get the available variables:


# 
# get_vars_for_dataset(dataset = "WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m")
# 


# # We can now download "bio01" and "bio02" for that dataset with:
# 
download_dataset(
  dataset = "WorldClim_2.1_HadGEM3-GC31-LL_ssp245_2.5m"
)


# 
# The datasets are averages over 20 year periods (2021-2040, 2041-2060, 2061-2080, 2081-2100). 
# In pastclim, the midpoints of the periods (2030, 2050, 2070, 2090) are used as the time stamps. All 4 periods are automatically downloaded for each combination of GCM model and SSP, and can be selected as usual by defining the time in region_slice.
# 
# future_slice <- region_slice(
#   time_ce = 2030,
#   dataset = "WorldClim_2.1_HadGEM3-GC31-LL_ssp245_2.5m"
# )


# 
# # Alternatively, it is possible to get the full time series of 4 slices with:
# 
# future_series <- region_series(
#   dataset = "WorldClim_2.1_HadGEM3-GC31-LL_ssp245_10m",
#   bio_variables = c("bio01", "bio02")
# )
# 



# read rasters ------------------------------------------------------------
library(terra)

setwd(file.path(mywd, "Env_var", "Future_climate", "worldclim_2.1"))
dir()

bio_ssp245_2030 <- rast("wc2.1_2.5m_bioc_HadGEM3-GC31-LL_ssp245_2021-2040.tif")
bio_ssp245_2030 <- terra::crop(bio_ssp245_2030, mystudyarea4326, mask=TRUE)
bio_ssp245_2030 <- project(bio_ssp245_2030[[-c(8,9)]], crs("EPSG:3035"), res=1000)
names(bio_ssp245_2030) <- gsub("wc2.1_30s_", "", names(my_bioclim))

bio_ssp245_2030 <- resample(bio_ssp245_2030, my_bioclim[[1]])

plot(bio_ssp245_2030[[1]])

plot(bio_ssp245_2030[[1]]-my_bioclim[[1]])




# save
dir.create(file.path(mywd, "Env_var", "Future_bioclim"))

setwd(file.path(mywd, "Env_var", "Future_bioclim"))

for (i in 1:nlyr(bio_ssp245_2030)) {
  
  layer_name <- names(bio_ssp245_2030)[i]
  
  writeRaster(
    bio_ssp245_2030[[i]],
    filename = paste0(layer_name, ".asc"),
    filetype = "AAIGrid",
    overwrite = TRUE
  )
}



# 
# 
# # # High resolution layers -------------------------------------------------------------------------
# library(terra)
# 
# grassland <- rast("G:/Copernicus/Grassland_2023_10m_Iberia.tif")
# 
# # Set 255 to NA
# NAflag(grassland) <- 255
# 
# grassland <- resample(grassland, my_bioclim[[1]])
# grassland <- terra::mask(grassland, my_bioclim[[1]])
# 
# setwd(file.path(mywd, "GIS"))
# writeRaster(grassland, filename="grassland.asc", overwrite=T, NAflag=-9999)
# 
# 
# forest <- rast("G:/Copernicus/Forest_Type_2021_10m_Iberia.tif")
# # Set 255 to NA
# NAflag(forest) <- 255
# 
# forest <- resample(forest, my_bioclim[[1]])
# forest <- terra::mask(forest, my_bioclim[[1]])
# 
# setwd(file.path(mywd, "GIS"))
# writeRaster(forest, filename="forest.asc", overwrite=T, NAflag=-9999)
# 
#